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ORIGINAL ARTICLE 

Influence of a priori Information, Designs, and 
Undetectable Data on Individual Parameters Estimation 
and Prediction of Hepatitis C Treatment Outcome 

THT Nguyen, 1 - 2 J Guedj, 1 - 2 J Yu, 3 M Levi 4 and F Mentre 1 - 5 

Hepatitis C viral kinetic analysis based on nonlinear mixed effect models can be used to individualize treatment. For that 
purpose, it is necessary to obtain precise estimation of individual parameters. Here, we evaluated by simulation the influence 
on Bayesian individual parameter estimation and outcome prediction of a priori information on population parameters, viral 
load sampling designs, and methods for handling data below detection limit (BDL). We found that a precise estimation of both 
individual parameters and treatment outcome could be obtained using as few as six measurements in the first month of therapy. 
This result remained valid even when incorrect a priori information on population parameters was set as long as the parameters 
were identifiable and BDL data were properly handled. However, setting wrong values for a priori population parameters could 
lead to severe estimation/prediction errors if BDL data were ignored and not properly accounted in the likelihood function. 
CPT: Pharmacometrics & Systems Pharmacology (2013) 2, e56; doi:1 0. 1 038/psp.201 3.31 ; published online 17 July 2013 



Chronic infection with hepatitis C virus (HCV) is a liver dis- 
ease that affects about 150 million people worldwide and is 
directly responsible of about 350,000 deaths every year. 1 The 
goal of anti-HCV treatment is to achieve a sustained virologic 
response (SVR), defined as undetectable serum HCV RNA 
24 weeks after treatment cessation. 2 HCV is classified into 
six major genotypes (GT) with HCV GT-2/3 being the sec- 
ond cause of chronic hepatitis C (after GT-1), accounting for 
-15-20% of infection in Western countries. 3 

Since 2001 , the combination of pegylated interferon (peg- 
IFN) and ribavirin (RBV) is the backbone of anti-HCV treat- 
ment with SVR rate of -50 and 80% in patients infected with 
GT-1 and GT-2/3, respectively. 4 - 6 In 201 1 , the approval of two 
protease inhibitors marked a new era of HCV therapy with a 
dramatic improvement in SVR rates in patients infected with 
HCV GT-1. 7-13 However, there is no clear evidence that Pis 
are beneficial in GT-2/3 patients. 14-17 Even though new treat- 
ment, such as nucleotide analogs may be effective against 
GT-2/3, 18 their cost and the fact that peg-IFN/RBV is already 
efficient makes bitherapy likely to remain essential in the 
treatment against HCV GT-2/3. 1517 Because peg-IFN/RBV 
therapy is associated with several significant side effects and 
high costs, 19 several efforts have been made to evaluate the 
possibility of treatment individualization. 20 

For that purpose, one can use viral kinetic models whose 
parameters have a high predictive value of treatment out- 
come. 2122 However, the use of these models is limited by the 
fact that frequent viral load data in the first weeks following 
treatment initiation are required to obtain precise estimation 
of the parameters. One way to improve the precision of indi- 
vidual parameter is to consider that the population param- 
eters are known a priori and to perform Bayesian estimation 
of individual parameters. Thus, this method combines a priori 



information of population parameters gathered from previous 
studies and individual viral load data prospectively obtained 
in a patient. This approach is similar to what is done in thera- 
peutic drug monitoring using population pharmacokinetic 
models. 2324 

However, the relevance of this approach is still contingent 
on the study design. 25 In practice, the difficulty to frequently 
assess viral load levels often gives predictions based on a lim- 
ited number of viral load data within each patient. We would 
therefore like to evaluate the quality of individual parameter 
estimation and treatment outcome prediction using a realistic 
design based on a small number of short-term observations. 

A common challenge in analyzing HCV kinetic data is the 
fact that a large proportion of viral load data are below the 
detection limit (BDL). Several studies have shown that naive 
approaches that omit or impute BDL data at an arbitrary value 
led to biased population parameter estimates, and this can 
be corrected by taking BDL data into account in the likelihood 
function. 26-29 However, these studies focused on the popula- 
tion parameters and did not evaluate whether and how these 
methods improved Bayesian individual parameter estimation. 

Here, our goal is to evaluate, by simulation and in the con- 
text of HCV GT-2/3, the influence of a priori population model, 
of viral load sampling designs and of methods for handling 
BDL data when estimating individual parameters and treat- 
ment response. 

RESULTS 

Description of the simulated data 

The percentages of BDL data were equal to 57.4, 27.8, 37.7, 

and 38.5 with designs D 0/l ,0^,0, , and D^ . „ , reSpeC- 
cs 24w' 4w' 4w_sparse' 4w_cha\\enge l ~ 

tively (see "Methods"). The SVR rate was equal to 78.4%, con- 
sistent with the rate of 80% reported in the literature for HCV 
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GT-2/3. 5 ' 6 Of note, 53.6% of the simulated patients had already 
eradicated the virus by week 12 (W12) of treatment. Non-SVR 
patients were classified as relapsers (undetectable HCV RNA 
at the end of treatment becomes detectable 24 weeks later), 
null-responders (HCV RNA declines from baseline <1 log 10 
at all time), or rebounders (HCV RNA declines from baseline 
>1 log 10 then rebounds under treatment), accounting for 7.6, 
2.3, and 1 1 .7% of our simulated population, respectively. Fur- 
thermore, two types of rebounders were considered: patients 
having a breakthrough (rebounds after a period being unde- 
tectable) or a partial virologic response (rebounds after a 
declining period but never being undetectable), accounting for 
3.3 and 8.4% of the simulated population, respectively. 

Individual parameter estimation 

Table 1 gives the mean relative bias (MRE), relative root 
mean square error (RRMSE), and shrinkage of the param- 
eters (5, 8, c, and e obtained for each design using true (M true ) 
or false population models (M 5e , M p ) and methods that omit 
(OMIT) or handle BDL data in the likelihood function (LBA). 

Even when the a priori information on population param- 
eters was correct (i.e., M true ), individual parameters were 
estimated with some bias and with the most informative 
design (D 24w ), RRMSE for the infectivity rate p was still 



equal to 55.4%. Moreover, p consistently had an extremely 
high shrinkage, >80% regardless of designs (Table 1 and 
Figure 1). This is due to the fact that this parameter can be 
precisely estimated only if there is a virologic rebound, a fea- 
ture that was observed in only 1 1.7% of patients. The three 
other parameters c, 8, and £ could be estimated with a shrink- 
age and a RRMSE <50% as long as viral load data in the 
first week of treatment were available. As expected, the preci- 
sion of parameter estimation deteriorates with less frequent 
sampling designs and the RRMSE for 8, c, and log 10 (1-£), 
increased by -30, 60, and 15%, respectively, when using the 
D 4w sparse as opposed to the D 4w design. Similarly, the shrink- 
age" for 8, c, and £ increased by -80, 90, and 300%, respec- 
tively, when using the D 4w sparse as opposed to the D 4w design. 
Surprisingly, a correct handling of BDL data did not critically 
improve the estimation of individual parameters when true 
a priori population parameters were set: for all parameters, 
only a small difference was found between LBA and OMIT 
methods using MRE (<6%) and RRMSE (<4%) as evaluation 
criteria, regardless of designs. 

Setting the fixed effect of 8 and £ at false values (M 5e ) 
naturally led to larger estimation errors, especially for 8 and 
c. Since in this model the a priori value for 8 was less than 
half of the true value (0.14 vs. 0.32 day -1 , respectively), 



Table 1 Mean relative bias, RRMSE, and shrinkage (in %) of /3, 8, c, efor different scenarios 3 



A priori 
information 


Design 


Method 




MRE (%) 






RRMSE (%) 






Shrinkage (%) 




P 


8 


c 


iog 10 d- 


-e) P 


8 


c 


"og 10 (i-£) 


P 


5 


c 


8 


M lrue 


D 24w 


LBA 


12.1 


4.9 


-3.5 


3.9 


55.2 


26.0 


27.9 


46.4 


89.8 


22.1 


27.5 


7.8 






OMIT 


13.2 


0.5 


-0.5 


1.9 


56.4 


25.7 


31.0 


47.6 


89.5 


28.8 


28.4 


23.3 




D 4w 


LBA 


12.4 


5.3 


-3.6 


4.3 


57.6 


27.9 


28.0 


41.2 


95.1 


26.0 


29.1 


8.9 






OMIT 


12.9 


1.0 


-0.8 


2.5 


58.3 


29.8 


31.0 


42.8 


94.7 


33.7 


29.9 


24.7 




^4w_sparse 


LBA 


12.2 


5.8 


6.0 


3.9 


57.5 


36.7 


45.4 


47.2 


96.3 


45.8 


55.6 


26.5 






OMIT 


12.4 


0.0 


8.9 


-0.1 


57.6 


35.0 


47.7 


50.5 


96.8 


51.5 


56.3 


48.7 




^4w_challenge 


LBA 


13.2 


3.9 


8.2 


4.9 


58.7 


41.6 


48.4 


65.2 


99.3 


61.8 


63.1 


67.2 






OMIT 


13.2 


3.5 


8.4 


4.4 


58.8 


41.6 


48.6 


65.3 


99.3 


60.0 


64.0 


68.2 




D 24w 


LBA 


12.8 


-7.2 


5.1 


2.3 


55.5 


23.7 


32.3 


29.6 


89.1 


27.1 


14.8 


1.7 






OMIT 


13.6 


-14.3 


11.8 


4.2 


56.5 


28.4 


41.6 


57.4 


89.1 


31.2 


12.2 


13.3 




D 4w 


LBA 


11.5 


-8.6 


6.0 


8.0 


56.4 


25.1 


33.0 


46.6 


95.6 


26.6 


16.2 


5.2 






OMIT 


11.9 


-15.7 


12.7 


5.7 


57.0 


29.6 


42.2 


49.4 


95.9 


33.6 


13.6 


16.7 




D„ 

4w_sparse 


LBA 


11.4 


-14.3 


27.6 


12.2 


56.7 


32.1 


63.0 


55.9 


97.6 


46.1 


54.8 


21.9 






OMIT 


11.9 


-24.6 


35.3 


-2.2 


57.3 


38.1 


70.2 


64.8 


97.8 


56.9 


53.5 


22.1 




^4w_challenge 


LBA 


11.5 


-20.7 


30.2 


-13.5 


57.6 


38.4 


65.1 


62.0 


99.7 


77.9 


54.4 


76.1 






OMIT 


12.6 


-35.1 


42.1 


-4.1 


58.4 


45.9 


76.1 


72.7 


99.8 


69.6 


61.5 


87.5 




D 24w 


LBA 


1028.3 


112.2 


5.1 


97.6 


1161.6 


123.5 


32.3 


228.3 


89.1 


27.1 


14.8 


1.7 






OMIT 


1036.8 


95.9 


11.8 


95.0 


1172.7 


111.1 


41.6 


228.3 


89.0 


31.2 


12.2 


13.3 




D 4w 


LBA 


1015.2 


109.0 


6.0 


99.3 


1155.6 


121.7 


33.1 


221.1 


95.6 


26.6 


16.2 


5.2 






OMIT 


1018.9 


92.6 


12.7 


97.0 


1161.3 


109.0 


42.2 


221.7 


95.9 


33.6 


13.6 


16.7 




^4w_sparse 


LBA 


1013.6 


95.9 


27.6 


103.8 


1156.0 


116.3 


63.1 


229.6 


97.6 


46.1 


54.8 


21.9 






OMIT 


1019.2 


72.2 


35.3 


89.1 


1163.1 


98.1 


70.2 


229.1 


97.8 


56.9 


53.5 


22.1 




^4w_challenge 


LBA 


1014.9 


81.4 


30.2 


77.9 


1161.4 


109.9 


65.1 


225.4 


99.7 


77.9 


54.4 


76.1 






OMIT 


1026.0 


48.4 


42.1 


52.1 


1173.9 


83.2 


76.1 


220.8 


99.8 


69.6 


61.5 


87.5 



a Obtained using different sets of a priori information, different designs and methods for handling BDL data. 
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Figure 1 Scatter plots of estimated values vs. simulated values for p, 8, c, s under five selected cases (top to bottom): M true _D 24w JLBA, 
M tr D. ^LBk, M, D. cnflrc _LBA, M, D. cnar OMIT, and M R D. __LBA. The dotted line is the identity line. Lower left corner graph 

true — 4w sparse — ' be — 4w_sparse — ' be — 4w_sparse — ' ji — 4w_sparse — J or- 

has a different scale on the y-axis compared to the other plots as in this scenario, p was fixed at a value that was 1 0 times higher than the 
true value to estimate individual parameters. 



the individual values for this parameter were consistently 
underestimated with an increasing relative bias, compared 
with that obtained with true information: for instance, using 
the design D 4wchajjen ^ relative bias going from 3.9% with 
true a priori population parameters up to 20.7% with false 
a priori information even when correctly handling BDL data. 
As the rate of the first phase of viral decline is proportional to 



cxe (see "Methods"), the estimates of these two parameters 
were correlated. This explains why setting a lower £ than the 
one used to generate the data (90% instead of 99.6%) led 
to an overestimation of c. This was particularly visible when 
considering designs without frequent samplings where c is 
poorly identifiable. For instance, the relative bias increased 
from 6% to 28% and RRMSE increased from 45% to 63% 
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when comparing the D 4w sparse design under the true and false 
model. However, the quality of parameter estimates remained 
reasonable (RRMSE < 50%) when using designs with fre- 
quent sampling in the first week, at least as long as BDL 
data were correctly handled. Indeed, and unlike what was 
seen when using the true population parameters, omission 
of BDL data dramatically deteriorated parameter precision. 
This was particularly noticeable for 8 where MRE doubled 
when omitting instead of correctly handling BDL data: for 
instance, using D 4w design, MRE increased from -8.6% with 
LBA method to -15.7% when BDL data were omitted. 

However, the result that setting false population parame- 
ters is not critical to achieve precise parameter estimates (at 
least as long as BDL data are correctly handled) becomes 
wrong if parameters are poorly identifiable. Indeed, extremely 
high errors were obtained for all parameters when setting /3to 
a value 10-fold larger than the true one (M p , Table 1). 

Prediction of treatment outcome 

We then investigated the impact of individual parameter esti- 
mation on treatment outcome prediction for the simulated 
patients. Results are presented in Figure 2. The influence of 
simulated individual parameters on treatment outcome was 
assessed in Supplementary Material 1 online. 

In case of correct a priori information, the prediction of indi- 
vidual treatment outcome was very good: even with the most 
challenging design where individual parameter estimates 
were biased, the misclassification rate was lower than 6%. 
The choice of designs did not have much impact on the pre- 
diction of treatment outcome although a slightly better pre- 
diction was obtained with rich design. With less information 
(by shortening or reducing sampling frequency), more false- 
positive (FP) responses were obtained as the prediction 
tended to be influenced by the response of the profile using 
the mean parameter values, i.e., SVR with this parameter 
regime (Supplementary Material 2 online) Of note, the rel- 
evance of this approach can be examined by comparing mis- 
classification rate using viral kinetic model with what would 
be obtained if no information was available. In that case, one 
would predict SVR in all patients, leading to an error rate of 
21 .6%, in which all are FP. 

In case of false a priori information about 8 and e, we 
found, as previously, that a correct handling of BDL data 
was necessary to obtain good prediction and omitting BDL 
data dramatically increased false-negative (FN) responses 
(Supplementary Material 2 online), especially when 
designs become less informative: for instance, in compari- 
son with a correct handling of BDL data, FN increased from 
1 .6% to 23.5% and 0.5% to 71 .5% when omitting BDL data 
in designs D 4w _ sparse and 

If one would predict SVR using rapid virologic response 
(RVR), defined as undetectable HCV RNA at W4 after treat- 
ment initiation, the misclassification rate would be still very 
good (6.9%). This, actually, means that the use of a viral 
kinetic model and data in the first month improved treatment 
outcome prediction by only 1 .4% and 1 .2% if a priori informa- 
tion was true or not, respectively. Thus the use of viral kinetic 
models does not bring, with this parameter regime, a critical 
benefit. This may be due to the fact that we evaluated the 
ability of treatment outcome prediction in a very favorable 



case where treatment outcome was highly correlated with the 
virologic response at W4. This is why we also compared the 
predictive ability of RVR and modeling in the case of shorter 
treatment duration of 1 2 weeks. In this more challenging case, 
the superiority of modeling approach was clear, with misclas- 
sification rates of only 14.6% and 19.2% when using correct 
or false a priori information, respectively, as compared with 
50% when using RVR (Supplementary Material 3 online). 

Whereas using M 5e did not much deteriorate the prediction 
quality, assuming false a priori information on p led to very 
high error rates (>20%) for all designs, even when BDL data 
were correctly handled, which is the consequence of a poor 
estimation of all individual parameters (see Supplementary 
Material 1 online). 

We illustrated the influence of different factors on indi- 
vidual estimation and outcome prediction using exam- 
ples of five typical patterns of virologic responses found 
in the simulated population (Figure 3). As expected, long 
design allowed obtaining individual fits close to the true, 
i.e., simulated profiles, especially when a virologic rebound 
occurred. With lower d and £, it is not surprising that M 5e 
predicted higher treatment duration to achieve SVR, com- 
pared to the true model. This explained high rate of FN 
obtained with this false information (Supplementary Mate- 
rial 2 online). Similarly, a larger number of treatment fail- 
ures were predicted with higher p, which is associated with 
more chances of virologic rebound. 

DISCUSSION 

Here we evaluated the role of a priori information, designs, 
and methods for handling BDL data in the prediction of 
individual virologic response and thus the ability of using 
viral kinetic models to individualize therapy in HCV GT-2/3 
infected patients. Unlike what is usually done, we focused 

75 -i 1 



72 - 




Figure 2 Misclassification rates (sum of false-positive (filled with 
pattern) and false-negative rates) of prediction of treatment outcome 
for a 24-week treatment in different scenario. 
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Figure 3 Examples of some individual fits for the different patterns of virologic responses obtained with different models, designs and BDL 
data handling methods:(a) SVR, (b) relapse, (c) null-response, (d) breakthrough, (e) partial response and rebound. The solid lines represent 



the predicted trajectories with the designs M D LBA (blue), M D 



LBA (magenta), 



-/> ^ 4w_sparse—'~ 

OMIT (green), and M^_D 4w spars ^LBA (red) and the true profile (black). Vertical line indicates the end of treatment (24 weeks). Observed and 
BDL data are presented as blue and red diamonds, respectively. Data that were not included in sparse designs are circled. 4-week designs do 
not contain measurements after W4. Simulation is stopped when the cure boundary is achieved, indicating that the patient has SVR. 



LBA (orange), 



^Se— ^4w_sp 



on individual parameters and assumed that the population 
parameters (i.e., the mean and the between-patient vari- 
ability) were known to estimate individual parameters using 
Bayesian approach. 

Using true a priori population parameters, precise esti- 
mation of individual parameters and good prediction of 
treatment outcome were obtained. Obviously, a major limi- 
tation of Bayesian approaches is that a priori values of the 
population parameters are required. This a priori informa- 
tion could be wrong, for instance if coming from a set of 



patients with different characteristics. This is why we also 
evaluated the robustness of this approach when the popu- 
lation parameters were different from the "true" parameters, 
i.e., those used for data simulation. Our results showed that 
setting wrong population parameters did not dramatically 
deteriorate individual parameters estimation as long as the 
parameters were identifiable (such as 5 and e). Importantly, 
satisfactory parameter estimation was obtained using a clin- 
ically relevant design with a limited number of data points 
(<6) and a last data point at W4. Therefore, the estimation 
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of viral kinetic parameters with this approach could allow 
for a rapid prediction of treatment outcome. In the case of 
p, a parameter that can be correctly estimated at individual 
level only in patients who have a virologic rebound (11.7% 
of our simulated population), assuming a higher value 
led to very large estimation errors in all individual param- 
eters. The combination of these largely biased parameters 
led to high prediction errors for treatment response (see 
Supplementary Material 1 online). 

BDL data is another factor that could lead to substan- 
tial bias. Interestingly, ignoring BDL data did not lead to 
significant bias as long as a priori population parameters 
were true. However, a correct handling of undetectable data 
became critical when population parameters were wrong, 
and we showed that this can be achieved by likelihood-based 
approaches that rigorously calculate the information con- 
tained in BDL data. 

Here we predicted the treatment outcome in a patient 
using the individual parameter estimates, which are the 
modes of a posteriori distribution, and therefore defining a 
yes/no response. In order to apply this approach to clinical 
practice, one would need to evaluate uncertainty of treatment 
outcome prediction by using the whole a posteriori distribu- 
tion of individual parameters, i.e., by incorporating uncer- 
tainty on individual parameter estimates. This could be done 
by reconstructing the a posteriori distribution, for instance 
using MCMC methods. 

The ultimate goal of viral kinetic modeling is to allow for 
individual parameter estimation and treatment tailoring. To 
achieve this goal, pharmacometricians will need reliable 
models and appropriate sampling designs that can provide 
precise parameter estimation. The information contained in 
a design can be a priori evaluated and optimized using the 
Fisher information matrix. This approach became a standard 
tool for pharmacokinetic/pharmacodynamic studies using 
NLMEM and was implemented in several software such as 
PFIM, PopED, POPDES, and POPT. 30 However, viral kinetic 
models have specific features that are not always properly 
addressed in these software such as the large proportion of 
BDL data 3132 and the need to achieve a good level of preci- 
sion of both population and individual parameters. 33 34 Clearly, 
optimal design at population and individual levels may dif- 
fer, and theoretical studies are needed to provide relevant 
designs at both levels. 

Here we studied the ability of modeling approach in the 
prediction of IFN-based treatment outcome in HCV GT2/3 
patients using a standard viral kinetic model. More sophisti- 
cated models are now being developed to explain complex 
viral kinetic patterns observed with new direct acting anti- 
viral agents. 35 Further studies will be necessary to suggest 
relevant designs to precisely estimate treatment outcome 
with this new generation of treatments and models. 

METHODS 

Viral kinetic model. In this study, we used a standard HCV 
model with liver regeneration 3637 given by the following set 
of equations: 



dt { T max j 

f t =pVT-8l (1) 

— = h-s)pi-cV 
dt v )H 

The model considers target cells T, infected cells I, and free 
virus V. Targets cells are produced at rate s, proliferate with 
rate r, die at rate c/and become infected at rate p. Infected 
hepatocytes die at rate 8 and release virions at rate p per 
cell. Virions are cleared from serum with rate c. Peg-IFN acts 
by blocking virion production, p, with effectiveness e varying 
between 0 and 1 . We assumed that the viral infection sys- 
tem before treatment initiation (t = 0) is at steady state. If 
only viral load data are measured, the parameters s, r, d, 
T max , p cannot be identified individually and therefore were 
usually set to values found in literature in parameter estima- 
tion. 22 ' 2538 In this study, we investigated only four identifiable 
parameters p, 8, cand e, which are also the most clinical rel- 
evant parameters. Indeed, in most patients treated with peg- 
IFN, viral load initially declines in a biphasic manner. This 
biphasic decline is captured by c, 8 and e. The first phase 
has a rate c x £ and leads to a viral load decline equal to 
log 10 (1-£), which reflects drug effectiveness. Subsequently, 
viral load declines with a slower but persistent rate approxi- 
mately equal to 8x e. If e is high enough, viral load will con- 
tinuously decrease during treatment. Otherwise, a virologic 
rebound may occur, 37 and the kinetics of viral rebound will 
depend on the value of p: higher values of p give earlier 
rebounds. 

Statistical model. Let / denote the /th individual (/= 1 AQ 

and /the /th measurement of an individual (j= 1 n. where 

n. is the number of observations for individual /). The statisti- 
cal model for the viral load jA.measured in individual / at time 
t.. is given by: 

y f = f(f fl 0,)xexp(e f ) (2) 

where f: nonlinear viral kinetic model, identical for all 
individuals, 
6:. individual parameters vector, 
e..: residual error supposed to follow normal distribu- 
tion N(0,cf). 

Because viral kinetics is simulated on long time period, it 
can happen that viral load takes values very close to 0. In 
that case, numerical errors may give a negative value for the 
viral load. This is why, unlike what is usually done, we evalu- 
ated the viral load in its natural scale and not in the log 10 scale. 

The individual parameters 0. can be decomposed into 
fixed effects /i = {p, 8, c, s} representing mean effects of 
the population and random effects r\ j specific for each indi- 
vidual. It is assumed that r\ j ~ A/(0, Q) with Q defined as the 
variance-covariance matrix in which, each diagonal ele- 
ments co q representing the variance of the qth component 
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of the random effect vector 77.. We assumed log-distribution 
for p, 8, c (qth parameter): 

e iq = ^ q exp(r liq ) 



and logit-distribution for e: 



(3) 



(4) 



' ^+(1-^)exp(-77 /e ) 

The population parameters vector y/\s composed of {u, Q, 0} 

Bayesian individual parameter estimation. In NLMEM, indi- 
vidual parameters are estimated using Bayesian approach. 
They are calculated as the Maximum a posteriori (MAP), i.e., 
the mode of the posterior distribution. Let y/ be the known 
population parameters vector and let Pfe I y,->v) be the con- 
ditional distribution of r\.. The MAP estimate of rj j is given by 



^ =argmax ri , (pfaly,^)) 

'p(y / h/^)xp(^/k) > 



=argmax 7ij 



p(y) 



(5) 



As /i is known, once 77, is estimated, 0, can be easily 
calculated. 

In presence of data below limit of detection (LOD), the dis- 
tribution p{y\r\ n \i/) is calculated by: 



p(y,ln,.v) = n-7*( l °9(yi). | 08('('J,.',)).^)i y| > t «, 



; n^(log(/-OD),log(^(77,,/ s )),cT 2 )l 



(6) 



y« <LOD 



Table 2 Parameter values of the different viral kinetic models in response to 
treatment (M true , M 5e , M ) 





Unit 


M true 




M p 


s 


cell • ml -1 • day -1 


60,000 


60,000 


60,000 


d 


day -1 


0.001 


0.001 


0.001 


8 


day -1 


0.32 


0.14 


0.32 


c 


day -1 


9 


9 


9 


P 


vrion • day -1 


50 


50 


50 


P 


ml • virion-Vday -1 


10- 7 


10- 7 


10- 6 


£ 




0.996 


0.900 


0.996 


r 


day -1 


0.006 


0.006 


0.006 


7 max 


cell • ml -1 


1.3x10 7 


1.3x10 7 


1.3x10 7 


co s 




0.5 


0.5 


0.5 


co c 




0.5 


0.5 


0.5 


% 




0.5 


0.5 


0.5 


co e 




2.5 


2.5 


2.5 


a 




0.461 


0.461 


0.461 



Four parameters that were simulated with random effects are /3, d, c, e. Pa- 
rameters that were modified to generate false models (M §e , M p ) are presented 
in bold. 



where (/)(x,m,v) is the probability density function, and 
0(x,m,v) is the cumulative density function of normal distri- 
bution with mean m, variance v, evaluated at x. 

Simulation study. This simulation study aimed at mimicking 
the viral response observed in HCV GT-2/3 infected patients 
during peg-IFN/RBV. The values for £, 8, and c in GT-2/3 
infected patients are based on literature. 39 For p, this param- 
eter has never been estimated in this population, and there- 
fore we use a value being in the range of possible values and 
close to those estimated in HCV GT-1 infected patients. 3640 
We used population parameters presented in the M true col- 
umn (Table 2) for data simulation. 

Four nested designs D 4w chaHenge cz D 4wsparse e D 4w e D 24w 
were considered. The richest design, D 24w , contains 12 sam- 
pling times until the end of the treatment (days 0, 1,4, 7, 
W2, 4, 6, 8, 1 2, 1 6, 20, 24). D 4w (days 0,1,4, 7, W2, 4) does 
not have any data after W4, D 4w sparse (days 0, 7, W2, 4) does 
not have any data after W4 and during the first week and 

D 4w challenge nas onlv one measurement at baseline and W4 
(day 6, W4). 

We simulated N = 1 ,000 vectors of random effects r\ j to 
calculate individual parameters 0. from the population param- 
eters of the true model. We then predicted the viral kinetics in 
these 1,000 in silico patients using the mathematical model 
(Eq.1). Residual errors were then simulated and added to 
viral load data for all sampling times of D 24w \o obtain dataset 
for D 24w design. As the four designs are nested, to generate 
datasets for other designs, we selected the simulated viral 
loads corresponding to the sampling times of the design of 
interest (Supplementary Material 4 online). As the same 
data of the same 1,000 simulated individuals were used 
to generate all the datasets, all the differences obtained in 
results were not due to Monte Carlo errors. 

We considered a limit of detection LOD = 45 lU/ml. Since 
successful therapy leads to SVR, i.e., the eradication of HCV, 
we hypothesized the existence of a threshold, called "cure 
boundary," under which infection is considered eradicated. 
Following what has been done previously, we assumed 
that this threshold corresponds to having <1 infected cell 
in the whole extracellular body fluid, assumed to be 15L. 22 
Although this cure boundary has not yet proven biologically 
correct, prediction of treatment outcome using this threshold 
has been shown to match SVR rate observed during peg- 
IFN based therapy. 22 If the course of infected cells crossed 
the cure boundary, we considered that SVR was achieved 
and no viral relapse could occur afterwards. All viral loads 
obtained after achieving the cure boundary are considered 
undetectable. 

To evaluate the impact of a priori information on individual 
parameter estimation and outcome prediction, we modified 
fixed effects of some parameters to generate two false mod- 
els (Table 2): for M 5e , we set 8 and e to the values obtained 
in GT-1 infected patients; for M p we set p to the upper limit of 
possible values found in literature. 40 

To estimate individual parameters for 1,000 simulated 
patients, individual random effects r\ j were first estimated and 
then used to calculate individual parameters. We supposed 
that population parameters were already known and fixed 
them at different sets of a priori values corresponding to M , 
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M 5e , M p to perform MAR BDL data were handled using two 
approaches: taking into account BDL data in the likelihood 
function (likelihood-based approach or LBA) or omitting all 
BDL observations in the dataset (OMIT). 

To evaluate the ability of estimation of the component q of 
an individual parameter, we used the following criteria: Mean 
Relative Error (MRE) and Relative Root Mean Square Error 
(RRMSE) defined by:t 



Study Highlights 



61,0-0,, 



MRE = — ^ — 



RRMSE q = 



1 N 



Oi,a-( 



(7) 



(8) 



For s, we evaluate the MRE and RRMSE of log 10 (1-e). 

We also used shrinkage (Shr) to evaluate the quality of indi- 
vidual parameter estimates (Eq. 9). A low shrinkage (<50%) 
indicates that individual parameter is correctly estimated. 



var 



Shr=1- 



(9) 



Once estimated, individual parameters were then used 
to simulate the individual predicted profiles of target cells, 
infected cells, and viral loads assuming a 24- or 12-week 
treatment. If infected cells were predicted to achieve the cure 
boundary during treatment, the predicted treatment outcome 
is having SVR and vice versa. 

By comparing the predicted treatment outcome with the 
true (simulated) response, we calculated the percentage of 
false-negative FN (patients having SVR predicted to remain 
infected), false-positive FP (uncured patients predicted to 
have SVR) and of misclassification MC (sum of FN and FP). 
Using these criteria, we compared the quality of treatment 
outcome prediction between different scenarios and with the 
case where no individual observation is available. To evaluate 
the gain of using modeling approach compared to empirical 
rules used in clinics, we also evaluated the predictive abil- 
ity for SVR of "rapid virologic response," RVR — defined as 
undetectable viral load at W4. We also considered a shorter 
treatment that was stopped at W12 in order to evaluate the 
ability of outcome prediction in a challenging case. 

All simulation was done using R 2.14.0. Parameter esti- 
mation (Supplementary Material 5 online) was performed 
using MONOLIX 4.1 .2. 
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WHAT IS THE CURRENT KNOWLEDGE ON THE 
TOPIC? 

S Viral kinetic models can be used for treatment 
outcome prediction in hepatitis C virus infected 
patient. Individual parameters can be estimated 
using Bayesian approach and a priori informa- 
tion on population model. 

WHAT QUESTION THIS STUDY ADDRESSED? 

y We evaluated the influence of a priori informa- 
tion, viral load sampling designs, and methods 
handling data below detection limit (BDL) on 
individual parameter estimation and treatment 
outcome prediction. 

WHAT THIS STUDY ADDS TO OUR KNOWLEDGE 

S Using correct a priori information, good indi- 
vidual parameter estimation and treatment out- 
come prediction were obtained with a "real-life" 
design containing only six measurements in 4 
weeks after treatment initiation. The results re- 
mained satisfactory if wrong a priori values on 
identifiable parameters were used unless BDL 
data were omitted. 

HOW THIS MIGHT CHANGE CLINICAL 
PHARMACOLOGY AND THERAPEUTICS 

1/ Bayesian individual estimation could provide 
rapid and reliable treatment outcome prediction 
to guide future therapy. 
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